clc;
clear;

example = 1;
switch example
    case 1
        node = [-1, -1;
            0, -1;
            1, -1;
            -1, 0;
            0, 0;
            1, 0;
            -1, 1;
            0, 1;
            1, 1];
        elem = [5, 4, 1;
            1, 2, 5;
            6, 5, 2;
            2, 3, 6;
            8, 7 ,4;
            4, 5, 8;
            9, 8, 5;
            5, 6, 9];
    case 2
        node = [0, -1;
            -0.707106781187, -0.707106781187;
            0.707106781187, -0.707106781187;
            0,-0.333333334358;
            -1,0;
            -0.707106781187, 0.707106781187;
            1, 0;
            0, 0.333333332309;
            0, 1;
            0.707106781187, 0.707106781187];
        elem = [8, 10,9;
            3, 7, 4;
            4, 8, 5;
            4, 2, 1;
            5, 2, 4;
            4, 1, 3;
            10, 8, 7;
            8, 6, 5;
            7, 8, 4;
            9, 6, 8];
end

Np = size(node,1);
Ne = size(elem,1);
Ned = Np+Ne-1;
elem_ed = zeros(Ne,3);
edge = zeros(Ned,2);
edge_el = edge;
for e = 1:Ne
    n1 = elem(e, 1);
    n2 = elem(e, 2);
    n3 = elem(e, 3);
    tempe = [n1, n2;
        n2, n3;
        n3, n1];
    if e == 1
        edge(1:3,:) = tempe;
        count = 3;
        elem_ed(e,:) = [1, 2, 3];
        edge_el(1:3,1) = 1;
    else
        count1 = 0;
        for m = 1:3
            flag = 1;
            for i = 1:count
                if tempe(m,:) == [edge(i,2), edge(i,1)]
                    flag = 0;
                    elem_ed(e,m) = -i;
                    edge_el(i,2) = e;
                    break;
                end
            end
            if flag == 1
                count1 = count1 + 1;
                new_idx = count + count1;
                edge(new_idx,:) = tempe(m,:);
                elem_ed(e,m) = new_idx;
                edge_el(new_idx,1) = e;
            end
        end
        count = count + count1;
    end
end


midpoint = 0.5*[(node(edge(:,1),1) + node(edge(:,2),1)), (node(edge(:,1),2) + node(edge(:,2),2))];

differences = node(edge(:,2),:) - node(edge(:,1),:);
point1_3 = node(edge(:,1),:) + differences/3;
point2_3 = node(edge(:,1),:) + 2*differences/3;


triplepoint = [point1_3; point2_3];
center = zeros(Ne,2);
for i = 1:Ne
    V = node(elem(i,1:3),:);
    center(i,:) = mean(V);
end

new_node = [node; triplepoint; center];
Nb = Np+Ned;

e1 = elem_ed;
e1(e1 > 0) = e1(e1 > 0) + Np;
e1(e1 < 0) = -e1(e1 < 0) + Np+Ned;
e2 = elem_ed;
e2(e2 > 0) = e2(e2 > 0) + Np+Ned;
e2(e2 < 0) = -e2(e2 < 0) + Np;
new_elem = [elem, e1, e2, (Np+Ned+Ned+1:Np+Ned+Ned+Ne)']

showmesh(new_node, new_elem, edge, "P1", Np);

bded = [edge_el(edge_el(:,2)==0,1), edge(edge_el(:,2)==0,:)]




%% showmesh
function showmesh(node, elem, edge, basis_type, Np)
Nb = size(node,1);
Ne = size(elem,1);
Ned = size(edge,1);
figure;
% mesh
h = patch("Faces", elem(:,1:3), "Vertices", node(1:Np,:));
set(h, "facecolor", "w", "edgecolor", "k");
axis("tight", "equal");
hold on;
% elem
center = zeros(Ne,2);
for i = 1:Ne
    V = node(elem(i,:),:);
    center(i,:) = mean(V);
end
text(center(:,1), center(:,2), int2str((1:Ne)'), "FontSize", 10, "FontWeight", "bold", "Color", "r", "HorizontalAlignment", "right");
% node
shift = 0.001;
range = (1:Np)';
plot(node(range,1), node(range,2), "k.", "MarkerSize", 10);
text(node(range,1)+shift, node(range,2)+shift, int2str(range), "FontSize", 10, "FontWeight", "bold", "Color", "k");
if strcmp(basis_type, "P1b") || strcmp(basis_type, "P2") || strcmp(basis_type, "P3")
    plot(node(Np+1:end,1), node(Np+1:end,2), "k.", "MarkerSize", 10);
    text(node(Np+1:end,1)+shift, node(Np+1:end,2)+shift, int2str((Np+1:Nb)'), "FontSize", 10, "FontWeight", "bold", "Color", "k");
end
% edge
midpoint = 0.5*[(node(edge(:,1),1) + node(edge(:,2),1)), (node(edge(:,1),2) + node(edge(:,2),2))];
text(midpoint(:,1), midpoint(:,2), int2str((1:Ned)'), "FontSize", 10, "FontWeight", "bold", "Color", "b", "HorizontalAlignment", "right");
hold off
end